perm filename CURFIT.F4[NET,GUE] blob
sn#041474 filedate 1973-05-14 generic text, type T, neo UTF8
subroutine curfit(x,y,sigmay,npts,nterms,mode,a,deltaa,
c sigmaa,flamda,yfit,chisqr)
implicit double precision (a-h,o-z)
dimension x(100),y(100),sigmay(100),a(10),deltaa(10),sigmaa(10),
c yfit(100)
dimension weight(100),alpha(10,10),beta(10),deriv(10),
c array(10,10),b(10)
11 nfree=npts-nterms
if(nfree) 13,13,20
13 chisqr=0.
goto 110
20 do 30 i=1,npts
21 if (mode) 22,27,29
22 if(y(i)) 25,27,23
23 weight(i)=1./y(i)
goto 30
25 weight(i)=1./(-y(i))
goto 30
27 weight(i)=1.
goto 30
29 weight(i)=1./sigmay(i)**2
30 continue
31 do 34 j=1,nterms
beta(j)=0.
do 34 k=1,j
34 alpha(j,k)=0.
41 do 50 i=1,npts
call fderiv(x,i,a,deltaa,nterms,deriv)
do 46 j=1,nterms
beta(j)=beta(j)+weight(i)*(y(i)-functn(x,i,a))*deriv(j)
do 46 k=1,j
46 alpha(j,k)=alpha(j,k)+weight(i)*deriv(j)*deriv(k)
50 continue
51 do 53 j=1,nterms
do 53 k=1,j
53 alpha(k,j)=alpha(j,k)
61 do 62 i=1,npts
62 yfit(i)=functn(x,i,a)
63 chisq1=fchisq(y,sigmay,npts,nfree,mode,yfit)
71 do 74 j=1,nterms
do 73 k=1,nterms
73 array(j,k)=alpha(j,k)/dsqrt(alpha(j,j)*alpha(k,k))
74 array(j,j)=1.+flamda
80 call matinv(array,nterms,det)
81 do 84 j=1,nterms
b(j)=a(j)
do 84 k=1,nterms
84 b(j)=b(j)+beta(k)*array(j,k)/dsqrt(alpha(j,j)*alpha(k,k))
91 do 92 i=1,npts
92 yfit(i)=functn(x,i,b)
93 chisqr=fchisq(y,sigmay,npts,nfree,mode,yfit)
if (chisq1-chisqr) 95,101,101
95 flamda=10.*flamda
goto 71
101 do 103 j=1,nterms
a(j)=b(j)
103 sigmaa(j)=dsqrt(array(j,j)/alpha(j,j))
flamda=flamda/10.
110 return
end